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Abstract. Time-dependent plasma codes are a natural extension of 
static nonequilibrium plasma codes. Comparing relevant timescales will 
determine whether or not time-dependent treatment is necessary. In this 
article I outline the ingredients for a time-dependent plasma code in a ho- 
mogeneous medium and discuss the computational method. In the second 
half of the article I describe recombination in the early Universe as a de- 
tailed example of a problem whose solution requires a time-dependent 
plasma code. 



1. Introduction 

Time-dependent plasma codes are required for any plasma where equilibrium 
is not maintained due to a time-dependent physical process. If the relevant 
physical process operates on a timescale t that is shorter than the timescale 
to reach equilibrium (of atomic ionization stages, for example) then the level 
populations must be followed in a time-dependent manner. For example, in an 
ionized gas with temperatures around 10 4 K the dominant timescale to reach 
equilibrium of the gas is the hydrogen recombination timescale, 

t rec = — - = 1.15^ 8 ng ^ours. (1) 

a(T e )n e 

Here a(T e ) is the recombination coefficient, T e is the electron temperature, n e 
is the electron density, ng is the electron density in units of 10 9 cm -3 , and £4 is 
the temperature in units of 10 4 K (Ferland 2000). There are many examples of 
timescales shorter than this such as shocks or rapidly varying radiation fields. 
Some astrophysical examples where time dependent plasma codes are useful 

are: 

• the recombination epoch 

• cooling of the first cosmological objects by H2 

• structure formation and evolution in the Universe 

• heating of the intergalactic medium 

• planetary nebulae and photodissociation regions 

• the interstellar medium. 

Time dependence always implies nonequilibrium level populations. In this 
article I only consider the optically thin case, where the level population's effects 
on the radiation field can be ignored or treated in a simple way. 
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In the second part of this article I discuss recombination in the early Uni- 
verse in detail as an example of an astrophysical problem whose solution re- 
quires a time-dependent plasma code. The time dependence is crucial in this 
case because the expansion of the Universe is much faster than the hydrogen 
recombination time. 



2. Time-dependent Equations 

2.1. Time-dependent Rate Equations 

The computational method in its simplest form involves solving one basic set of 
equations 

d n i,k 

^ -^populate ^depopulate: 

constrained by particle conservation 

Y j n i , k {t) = N k (3) 

i 

and charge conservation 

Y / J2 in iM t )- n S) = °- ( 4 ) 

k i 

Here represents the number density of atoms of species k in an ionization 
stage i, n e is the number density of electrons, and N k is the total number density 
of all ions of a species k. Note that there is one redundant equation for each k 
and one redundant equation for each i of a given k. 

The population and depopulation terms in equation [2] for a given species k 
and ionization stage i can be described as follows: 

Rpopuiate = ^j-i,fc(r^i i fc+r^i i fc+r^ ljfc )+rei + i i fc(r^ c ljfe +r^ 1)fc -i-r^ 1)fc ), (5) 

^depopulate = n k (Tf k + Tf k + T^ k + T™ k + T^ k rec + vf£ ec ). (6) 
Here V are the rate coefficients: 

Photoionization : Tf k = I —Oi^(i>)J{y,t)dv (7) 



Recombination : Y\^ k = n e aj_i ) fc(T e ) (8) 

C/C,rec _ f C/( 
i,k — n eJi,k 



Collisional ionization/recombination : Ff[. C ' rec = n e ffl C ' rec {T e ) (9) 



Charge — exchange ionization/recombination : 

r i = n H+ ^ k (T e ) + n He+1 ^{T e ) (10) 
rjT = nnlUT e ) + n He ^(T e ). (11) 
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Here J(v, t) is the radiation field, cr{v) is the photoionization cross section, a(T e ) 
is the recombination coefficient, and h is Planck's constant. The equations are 
coupled among different species by the charge exchange reactions. Not shown 
in the above set of equations are that the excited states j can be followed in the 
same framework, rather than only considering the ground state of an ionization 
stage. This can be time consuming because it involves j more equations for each 
i and k. One simplification is to only use the lower levels in the photoionization 
equations and to use a pre-computed recombination coefficient that takes into 
account recombination to the upper energy levels of the atom. Another simpli- 
fication is to construct several averaged atomic energy levels ("super levels"; see 
Lucy 2001) in lieu of hundreds of excited upper levels. This works because the 
individual excited states (above the first several n levels) usually don't domi- 
nate most rates but their collective effects are important. The rate equations 
can also be extended to include molecules, in which case the terms R p0 puiate and 
^depopulate involve reaction rates between different molecular species (see e.g., 
Stancil, Lepp, & Dalgarno 1998). 

When nonequilibrium populations are involved one must take into account 
their effects on the electron temperature. If the deviations from equilibrium are 
strong, the temperature may depart from its equilibrium value. Furthermore the 
time dependence is key. The radiation field is also affected by the nonequilibrium 
values although in plasma codes the radiation field is not solved for explicitly. 

2.2. Time-dependent Temperature Equation 

In a static situation the temperature is determined by the equilibrium between 
heating and cooling processes. In the time-dependent case the heating and 
cooling processes might not have had time to reach equilibrium and so they 
must be solved for together with the time-dependent rate equations. In general 
one can test whether or not a given heating and cooling process has reached 
equilibrium by comparing its timescale to the time-dependent process (shocks, 



variation of radiation field, wind velocity, etc.); see § 4.1 j for an example. 
The time-dependent kinetic energy equation is 

jE{t) = G(t) - L(t), (12) 
where G(t) are the energy sources and L{t) are the energy sinks: 

/°° Air 
—a i>k {v)J{y, t)h{v - v Q )dv, +G b \t) (13) 

= E E n i+ i, fc (t)n 6 (t)a i _a, fc (r e ) + L bb (t) + L^(t). (14) 

k i 

In this example the heating processes considered are photoionization heating 
and collisional line heating G bb (t) (the first and second terms in equation 13 



respectively) while the cooling processes are recombination cooling, collisional 
line cooling L bb (t), and bremsstrahlung cooling (t) (the first, second, and 
third terms in equation [l4| respectively) . For explicit descriptions of these terms 
see, e.g. Osterbrock (1989). While these processes tend to be the most important 
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heating and cooling mechanisms, their relative importance and the importance 
of other processes depends on the specific situation (i.e. T, density, etc.) In 
general it pays to be clever by checking which heating and cooling rates are 
important in advance. For example, heating and cooling processes due to less 
abundant elements can usually be ignored. The time dependence means that 
the temperature (=2E/3N to tkB) at a previous time, together with the equations, 
determines the temperature at the present time. 

2.3. The Radiation Field 

In time-dependent plasma codes the radiation field is a key ingredient because it 
governs the level populations (ionization stages and atomic levels) and temper- 
ature. It is not the goal of plasma codes to solve for the radiation field, rather 
it is considered a given parameter in an optically thin situation. The detailed 
solution of the radiative transfer problem is a much more complex, difficult, and 
different class of problem. Still, small changes to the radiation field can be con- 
sidered. For example, in a plasma with a central radiation source (e.g. an HII 
region) the radiation field can be written 



where the first term on the right side is the radiation field from the central source 
modified by absorption, and the second term is the diffuse radiation. Here R is 
the radius of the central star and r is the radial distance from the star center. 
Note that in the optically thin case the modification to the radiation field by 
exponential attenuation can easily be computed. Here the optical depth r is 
simply 



Jn k i 

Usually the excited levels j above the ground state do not contribute much to 
the optical depth. 

3. Computational Method 

3.1. Explicit vs. Implicit Methods 

The time-dependent rate equations and the time-dependent temperature equa- 
tion form a set of coupled first-order nonlinear differential equations. The equa- 
tions are first order with respect to the independent variable, time. The equa- 
tions are nonlinear because of recombination and collisional terms that involve 
the interaction of two species, and n e . The time-dependent plasma calcu- 
lation is an initial value problem where the initial values need to be specified 
and the goal is to determine the number densities and temperature at a later 
time. There are two general classes of numerical methods used to solve initial 
value problems of first-order ordinary differential equations, explicit methods 
and implicit methods. 

Consider the rate equations to be of the following form 




(15) 




(16) 



-t: 1 = fi(t,n 1 ,...,n N ,n e ,T) 



(17) 
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where is the known equation fa = Rp_ opu i a te ~ Rdepopuiate described in §ETl] In 



the limit of small timesteps equation 17 can be rewritten as 



Arn = f i (t,n 1 ,...,n N ,n e ,T)At. (18) 

The new value n new can be determined with the derivative fi to propagate the 
solution, 

n new = n i d + Am^d. (19) 

This is known as the explicit method and is also called Euler's scheme. More 
sophisticated versions of it that use derivatives at other positions, for example, 
are more useful. For example a popular algorithm is the Runge-Kutta method. 

Explicit methods can be unstable for stiff equations — those equations where 
different no- are changing on very different timescales. Stiff equations are com- 
mon in a set of rate equations when some rates (e.g. a chemical reaction rate or 
an ionization rate) are very different from other rates that control relevant num- 
ber densities. Taking a small enough timestep to satisfy the shortest timescale 
is almost always impractically time consuming. Taking a reasonable timestep 
can have disastrous consequences (such as values approaching infinity or oscil- 
lating values) for species that evolve on the timescale(s) much shorter than this 
reasonable timestep. The following typical example illustrates this nicely. 

Consider the equation 

An,, = -RmAt. (20) 

The new value of is 

n n ew = n id + An ifi id (21) 
n n ew = n i d (l - RAt). (22) 

In the limit of a large timestep, At, |n new) | — > oo, an obviously unphysical value 
and one that certainly does not satisfy the constraint of equation |3|. 

Implicit methods are commonly recommended to deal with stiff equations. 



The form of the equations is the same as equation |18| but the derivative / is 
evaluated at the new time, 

nnew — n Q id ~\~ An.j new . (23) 

For the above example, n new becomes 

n new = n old /(l + RAt). (24) 

In this case, for a large timestep n new — ► 0, the correct solution to equation |2(]. 
Implicit methods for a set of nonlinear equations (as opposed to a single linear 



equation like equation 2C) often involve the Jacobian which must be used in an 



iterative approach to solve for the set of values n new ^. But for the rate equations 



described in equations g through |10|, the Jacobian can be computed analytically 
in a straightforward way. 

There are many numerical methods that implement the forward or backward 
Euler scheme into more sophisticated algorithms. Many textbooks deal with stiff 
equations (e.g. Press et al. 1992; Lambert 1973) and have detailed discussions 
about error analysis, adaptive stepsizes, and more. 
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3.2. Numerical Considerations 

The time-dependent rate equations and the particle and charge conservation 
equations mean there is one redundant equation for each k and one redundant 
equation for each i for a given species k. A practical approach is to solve for 
the time-dependent rate equations (equation ^) and use the particle and charge 
conservation equations (equations |3| and as a check to test the accuracy of 
the solution. 

The input parameters for the set of rate equations and the temperature 
equation are abundances, accurate atomic data (e.g. atomic energy levels, os- 
cillator strengths, absorption cross sections, rate coefficients, etc.), and a good 
starting solution. A good starting solution is usually the static solution, which 
is the same set of equations (equation ^ with dni jdt = 0) . If there are too many 
equations for an algebraic solution, the starting solution is best obtained by a 
Newton-Raphson type scheme. The output is the population number densities, 
ni,k(t), and the electron temperature. 

Algorithms with adaptive stepsize control and that monitor internal errors 
are extremely useful. The errors are used to ensure the solution reaches a speci- 
fied accuracy. When debugging it is useful to follow the errors (i.e. write them to 
a file). If different n^fc have the same error this indicates the problem is related 
to equations that involve those two number densities (e.g. the derivatives or 
terms in Jacobian). In general implicit methods with adaptive stepsize control 
with a specified accuracy should suffice for a time-dependent plasma code. 

4. Recombination in the Early Universe 
4.1. Introduction 

Three hundred thousand years after the Big Bang the Universe became cool 
enough for the ions and electrons to form neutral atoms. This was the recom- 
bination epoch. The Universe expanded and cooled faster than recombination 
could be completed and a small but non-zero fraction of electrons and protons 
remained. This is referred to the ionized fraction at freezeout. Recombination 
in the early Universe was first studied by Peebles (1968), Zel'dovich et al. (1969) 
and others, and the basic framework has remained unchanged since that time. 
For a complete set of references to both the early papers and more recent pa- 
pers see Seager et al. (2000). The recombination epoch is a good example of 
a problem that must be solved by a time-dependent plasma code because the 
expansion timescale of the Universe is shorter than the recombination timescale. 

Figure 1 shows the difference between a time-dependent calculation (de- 
scribed later in § |4.3,| ) and a time-independent calculation from the Saha equi- 
librium equation. Here x e = n e /riH where n e is the number density of elec- 
trons and rin is the number of total hydrogen nuclei. The ionization fraction is 
used instead of number density because it gives a clearer picture of recombina- 
tion; expansion of the Universe changes the volume and hence number densities 
and the total number density depends on the cosmological model. Similarly, 
redshift (z) is used instead of time because it is independent of the cosmo- 
logical model. Redshift is the Doppler shift of light from a receding source: 
z = AX/Xq = v/c = H^r/c (for low redshifts), where v is the recessional veloc- 
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Figure 1. The ionization history x e = n e /ri}j as a function of redshift 
and time. The time on the upper x axis is approximate. The solid line 
is x e from the time-dependent calculation and the dashed line is x e 
from the time-independent Saha equilibrium equation. The dotted line 
is the ratio of neutral H to total H nuclei. The cosmological parameters 
used in this calculation are Q m = 0.25, A = 0.75, and f^/i 2 = 0.02. 



ity, r is the distance, and Hq is Hubble's constant. Redshift and time are related 
by dz/dt = —(1 + z)H(z) where H(z) is the Hubble parameter, which depends 
on several cosmological parameters. For a flat £l m = 1 universe, this relation is 

_ 2 1 
^ZH^l + zW 1, 

but is more complicated for other cosmologies (see, e.g., Peebles 1993). Basically 
redshift is a convenient parameter because it can be measured directly and the 
ignorance of cosmological parameters is hidden in usage of redshift instead of 
time. The radiation temperature from adiabatic cooling is Tr = Tb(l + z), 
where To = 2.728 (Fixsen et al. 1996). To illustrate the density-temperature 
parameter space of recombination, Figure 2 shows the electron density as a 
function of redshift and of radiation temperature for the same ionization history 
and cosmological model shown in Figure 1. 

A recent motiviation to revisit the early Universe recombination calculation 
is that the ionization fraction history (shown in Figure 1) is a basic determinant 
of the Cosmic Microwave Background (CMB) spatial anisotropics which will be 
measured to the 1% level with the MAPQ satellite after its launch in June 2001 
and later with the PlanckQ satellite. Figure 3 shows the angular power spectrum 



1 ht t p : / / map . gsf c . nasa. gov / 

2 http: //astro. estec.esa.nl/SA-general/Projects/Planck/ 



(25) 
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Figure 2. The same ionization history as shown in Figure 1 but with 
the number density function of redshift and of temperature 

(Tr = 2.728(1 + z), and T M = T R until z ^ 100). The solid and dashed 
lines correspond to n e from the time-dependent and time-independent 
calculations. The dotted line is the number density of neutral H atoms. 
The number densities decline throughout recombination because of the 
expansion of the Universe. 
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Figure 3. Angular power spectrum of CMB anistropies from COBE 
and recent experiments. Courtesy of Ned Wright (Wright 2001). 
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of CMB anisotropies, where the C;s are squares of the amplitudes in a spherical 
harmonic decomposition of anisotropies on the sky. They represent the power 
and angular scale of the CMB anisotropies by describing the rms tempeartures 
at fixed angular separations averaged over the whole sky (see e.g., White, Scott, 
& Silk 1994). These temperature differences of the CMB at fixed angular scales 
on the sky correspond to the "seeds" from which galaxies and other structures 
grew. Cosmological parameters (such as Qq, ^s, h, etc.) can be determined 
from fits (notably the relative heights and positions of the peaks) to the power 
spectrum. 

It is instructive to ask: when did the Universe become neutral? As a simple 
guess one might think that recombination occurred when the CMB radiation 
peak had just cooled to B\ = 13.6 eV, the binding energy of the ground state 
of H. This corresponds to a temperature of 56,000 K, which translates to a 
redshift of approximately 20,000, and a time after the Big Bang of roughly 3300 
years. A better guess would consider that there are many, many more photons 
than baryons (the ratio is ~10 9 ) in the Universe. This gives a temperature 
of ~7000 K or a redshift of ~2500 which corresponds to a Universe age of 
approximately 75,000 years. However, due to the short collisional timescale 
between H atoms (see Figure 4) the atomic structure of H is important, and the 
time of recombination is controlled by the plasma temperature. Therefore the 
most reasonable estimate comes from the Saha equilibrium equation, 



at the point where some fraction — say 99% — of the atoms have become neutral. 
The temperature at which this occurs is about 3000 K, which corresponds to a 
redshift of approximately 1000 and a time of approximately 300,000 years after 
the Big Bang. In this equation the UjS are partition functions, m e is the electron 
mass, ks is Boltzmann's constant, and the other constants are as previously 
defined. 

It is also useful to consider different timescales that are relevant during 
the recombination epoch, shown in Figure 4. In general the relevant timescales 
determine whether or not a time-dependent calculation is even necessary. The 
most important point in this example is that the expansion timescale (tn) be- 
comes shorter than the recombination timescale early on, meaning that a time- 
dependent treatment is crucial. Also, the recombination timescale is shorter 
than the Saha equilibrium timescale (tsaha), meaning that Saha equilibrium is 
not valid. The hydrogen atom collisional timescale (t co u), the Coulomb timescale 
(tcoul), and the proton collision timescale (= (m p /m e ) l / 2 tc ui) are all much 
shorter than the expansion timescale so matter (electrons, protons, H atoms) 
has had plenty of time to relax to a state of thermodynamic equilibrium. Thus 
we have Tm = T e . The Compton scattering timescale (tcomp) is faster than 
the expansion timescale, meaning that the photons and electrons remain at the 
same temperature throughout recombination. At z ~ 100 when tcomp ^ tn 
the matter (electrons) and radiation (photon) decouples. For the equations that 
describe these timescales see Scott (1988). From this example we see that the 
only equations that need to be followed in a time-dependent manner are rate 
equations (to follow the time dependence of recombination) and the matter tem- 
perature. 



n i+ in e _ 2u i+ i (2TTm e k B T M ) 3/2 
m Ui h 3 



-B 1 /k B T M 



(26) 
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Figure 4. Important timescales during the recombination epoch. See 
text for discussion. 

4.2. The Recombination Process 

Hydrogen recombination involves protons capturing electrons and electrons cas- 
cading down to the ground state. There are many recombinations to and from 
each energy level and many bound-bound transitions among the energy levels 
before there is one net recombination to the ground state. Many people op- 
pose the use of the term "recombination" to describe the very first time in our 
Universe's history that protons and electrons combined to form neutral atoms. 
There were many ionizations and recombinations for every net recombination, 
however, so the term recombination is still appropriate. 

The recombination process was not instantaneous (it was essentially case 
B but cf., Seager et al. 2000) because of the strong but cool CMB blackbody 
radiation field. The electrons, captured into different atomic energy levels, could 
not cascade instantaneously down to the ground state. The electrons were im- 
peded because of fast reionizations out of excited states that were due to the 
huge reservoir of low-energy photons and because of the high optical depth of 
the Lyman lines and the continuum transitions to the ground state. Any Lyman 
line or continuum transition to the ground state emitted a photon with energy in 
which there were few blackbody photons, and this immediately photoexcited or 
photoionized a neighboring atom in the ground state. Figure 5 shows a black- 
body radiation field with the energy levels of the 13.6 eV transition and the 
Lyman a transition where there are few blackbody photons. Atoms reached the 
ground state either by the 2s-ls two-photon process or through the cosmological 
redshifting of the Lya line photons. The cosmological redshifting occurs because 
as space expands the frequency of the photons changes — possibly enough to be 
redshifted out of interaction frequency with the Lyman a line. Because both 
of the rates from n=2 to the ground state (the 2s-ls two-photon process and 
the cosmological redshifting) were much slower than the net recombination rate 
to n=2, a "bottleneck" occurred that slowed down the entire recombination 
process. 
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0.01 0.10 1.00 10.00 100.00 1000.00 10000.00 

A (/xm) 

Figure 5. The solid line is the blackbody radiation intensity at 
z = 1200, which corresponds to Tr = 3274 K. The dashed line is 
the blackbody intensity at z = 500 or Tr = 1364. The symbols show 
the wavelengths (i.e. energies) of different H atom transitions. Here 
"cont to n=l" is the 13.6 eV continuum to ground state transition (the 
binding energy of the ground state H atom) . At the wavelengths of the 
continuum to n=l transition and the Lyman alpha transition — which 
are in the high-energy tail of the blackbody radiation field — there are 
very few blackbody photons. For the n=60 to n=30 transition (as 
an example) there are always plenty of blackbody photons throughout 
recombination and so the net rate of this transition is always the equi- 
librium one — zero. Other high-energy transitions such as n=70 to n=2 
are at a wavelength where there are many blackbody photons early on 
during the recombination epoch (here e.g. at z = 1200) but not at later 
times (here e.g. at z = 500). Thus these types of high energy tran- 
sitions go out of equilibrium at some point during recombination and 
cause a faster recombination (partly due to a faster cascade) compared 
to the equilibrium scenario. 
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4.3. Recombination Equations 

Recombination in the early Universe can be computed using the formalism de- 
scribed in §|2] 

Time- dependent rate equations We are interested in the time evolution of the 
H n-level populations, in addition to protons (ionized H) and electrons. We have 
one time-dependent rate equation for each atomic energy level of H, 



(1+*) 



dnj (z) 1 



dz H(z) 



N 

)P cj - nj {z)P jc ) + $>i^ 

3=1 



+3^-0), (27) 



and also equations for the proton number density (n p ) and electron number 
density (but if only H is being considered n p = n e so only one equation need be 
followed). Here ARji is the net bound-bound rate between bound states i and j, 
and the Pj c are the rate coefficients between bound levels j and the continuum 
c: Pij = Rij + n e Cij and Pj c = Rj c + n c Cj c , where R refers to radiative rates 
and C to collisional rates. Here the ns are physical (as opposed to comoving) 
number densities: rij refers to the number density of the jth excited atomic 
state, n e to the number density of electrons, and n c to the number density of 
a continuum particle (i.e. proton in this case). For convenience we use redshift 



z instead of time (see equation 25 and the accompanying discussion). H(z) is 
the Hubble factor and the extra term 3nj(z) comes from the fact that space 
is expanding. For a 300-level H atom — where the number of n levels needed 
is determined by a thermal broadening cutoff — there are 300 such equations 
and together they involve hundreds of bound-free transitions and thousands 
of bound-bound transitions. The number of n levels that should be considered 
depends on the thermal broadening cutoff. We do not consider individual I states 
(with the exception of 2s and 2p), but assume the I sublevels have populations 
proportional to (21 + 1). The £ sublevels only deviate from this distribution in 
extreme nonequilibrium conditions (such as planetary nebulae). Dell' Antonio & 
Rybicki (1993) looked for such i level deviations for n < 10 and found none. For 
n > 10, the £ states are even less likely to differ from an equilibrium distribution, 
because the energy gaps between the I sublevels are increasingly smaller as n 
increases. However the ^-level time evolution could easily be included in the 
plasma code but would be consume an unreasonable amout of computer time. 
The photoionization equation, 

f°° Air 

/ —aj(y)B(y,T R )dv, (28) 
and the recombination equation, 

( h2 V /2 a, R l/kB T M r ^ ■ ( 2hu 



2Trm e k B T M I 2g k 



Ej/kBTM r Vv^^ + B(y ' tr) ) e ~ kv/kBTMdv 



(29) 

are the familiar ones, but note the separate consideration of Tm and Tr. Here 
c is the speed of light, g is the statistical weight, E is the ionization energy 
and the other constants and variables have their usual meanings. Here we have 
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replaced the usual J(y,t) = B(v, Tr) (see the below discussion on the radiation 
field) because Tr is a function of z and hence t. Collisional rates are much 
smaller than radiative rates (due to the large reservoir of photons), so while 
they are usually important for a plasma of this temperature and density, they 
can actually be ignored in the early Universe recombination calculation. 

Every time a situation involves nonequilibrium level populations or ioniza- 
tion stage populations — and time-dependent cases are no exception — one must 
consider the effects of nonequilibrium populations on the temperature and on 
the radiation field. 

Time- dependent temperature equation The radiation temperature during re- 
combination is determined from adiabatic cooling of radiation, Tr = To(l + z), 
where Tq = 2.728 is today's temperature of the CMB determined from COBE 
measurements (Fixsen et al. 1996). The time-dependent matter temperature 
equation for recombination in the early Universe is 

{1 + z) dTM = SarU {T _ Tr) + 2T ^ (30) 

dz 3H(z)m e c n e + Nr- + Nu e 

where ax is the Thomson scattering cross section and U is the energy density. 
The first term on the right side is Compton cooling and the second term comes 
from adiabatic cooling of an ideal gas. Other heating and cooling terms such 



as ionization heating, bremsstrahlung cooling, and others mentioned in § 2.2 
are negligible and need not be considered. In general one can determine in 
advance of a calculation whether or not different heating or cooling terms are 
important from comparing timescales. For example in this case the Compton 
cooling timescale becomes longer than the expansion timescale at z ~ 100 (see 
Figure 4) and including Compton cooling in the time-dependent matter temper- 
ature equation keeps track of when Tm and Tr differ. Judging from Figure 4 it 
is not necessary to include the evolution of the matter temperature for a first 
order calculation, since it has little effect at z ;> 100. Nonetheless, the Tm and 
Tr difference is significant enough to affect the ionization fraction at freezeout 
by a few percent at low z. 

The Radiation Field The time-dependent radiation field is 

(1 + z) d J^A = 3J(l/j z)--^- \j(u, z) - k{u, z)J{v, z)] , (31) 

where z) is the emission coefficient and k(u, z) is the absorption coefficient. 
It would be extremely difficult and time consuming to solve this time-dependent 
radiative transfer equation at each redshift for many frequencies. During the 
recombination epoch the background radiation field was very smooth; it is gen- 
erally a blackbody radiation field and today the CMB is a blackbody as mea- 
sured by COBE. Furthermore the blackbody thermal spectrum is preserved in 
the expansion of the Universe; the time-dependent radiation field is simply a 
blackbody radiation field determined by Tr = 2.728(1 + z). This means that 



instead of solving equation 31 we can just use the blackbody intensity. However, 
the extreme trapping of Lyman line photons means there are significant distor- 
tions to the blackbody radiation field. In this case we can treat these distortions 



14 



Seager 



by the Sobolev escape probability — this is a solution of the radiative transfer 
equation in the presence of moving media, which in this case is the expanding 
Universe. This approach was first used by Dell' Antonio and Rybicki (1993) and 
makes the recombination in the early Universe problem tractable. The term 
ARij will depend on the escape probability from Sobolev theory for transitions 
where it is appropriate (the Lyman lines). 

Numerical Considerations A natural starting solution to the recombination in 
the early Universe problem is at high redshift (i.e. early times and high T) 
where all the hydrogen is ionized and the starting solution is n e = n p = Nh and 
rij = 0. Because not exactly all of the hydrogen is ionized at z ~ 2000 it is better 
to instead use the solution of the static rate equations using a Newton- Raphson 
scheme with the input to that as n e = n p = Nh and rij = 0. 

When using an integration algorithm with a specified accuracy it helps to 
not include the errors of very small level populations in the consideration of the 
size of the next timestep. I use the requirement to ignore the errors if nj/Nn ^ 
10~ 13 because populations that small are not relevant to the calculation anyway, 
but they are still important enough at other times to keep their number densities 
in the system of rate equations. Otherwise the accuracy requirements are too 
stringent and the integration will approach a very tiny stepsize and a very long 
computational time. 

4.4. The "Standard Recombination Calculation" 

The "stand ard" methodology forgoes the detailed numerical computation de- 
scribed in § |4.3, and considers an "effective three-level atom" with a ground 



state, first excited state (n=2), and continuum, with the n > 2 states repre- 
sented by a recombination coefficient. A single ordinary differential equation 
can then be derived from the rate equations shown in equation ^7| to describe 
the ionization fraction (see Peebles 1968, 1993; Seager et al. 1999). Many as- 
sumptions go into this derivation, including the following: that H excited states 
are in equilibrium with the radiation; that stimulated de-excitation is negligible 
for the Lya transition; that a simple recombination coefficient can be used; that 
every net recombination results in a ground-state atom, so that the ground-state 
number density n\ = Nh — n p ; that the Lya redshifting can be dealt with using 
a simple escape probability; that collisional processes are negligible; and that 
He can be ignored. Only the assumption that the upper levels of the H atom 
are in equilibrium with the radiation is not valid. This is described in the sub- 
section below. See Peebles (1968) or Seager et al. (1999) for the single ordinary 
differential equation that describes the ionization fraction history. 

4.5. The Exact Recombination History 

Figure 6 shows the recombination history for the "modern" detailed numeri- 
cal computation compared to the standard, single ordinary differential equation 
calculation. The modern, more detailed calculation has a faster recombination 
rate which results in a 10% smaller ionization fraction at freezeout compared to 
the standard calculation. This in turn has an effect on the CMB power spec- 
trum of anisotropies of a few percent. Note that the cosmological parameters 
(and other second order effects) determine the shape of the power spectrum, 
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Figure 6. Ionization histories from the standard recombination cal- 
culation compared to the more detailed "modern" calculation described 
in this article. 

and to derive them accurately the recombination history must be known to high 
accuracy. In the standard case that assumes equilibrium among the excited 
states n > 2 the net bound-bound rates are by definition zero, and this is an 
implicit assumption in deriving the single ordinary differential equation used in 
the standard calculation. We find that at z <; 1000, the net bound-bound rates 
become different from zero because, at low temperatures, the cool blackbody ra- 
diation field means that there are few photons for photoexcitation of high-energy 
transtions (e.g. 70-10, 50-4, etc.). This is shown in Figure 5 where the vertical 
lines show the wavelengths of various H atomic transitions, the solid line is the 
blackbody intensity near the beginning of recombination, and the dashed line 
is the blackbody intensity towards the end of recombination. Because there are 
few available photons for the high-energy transitions, spontaneous de-excitation 
dominates those bound-bound transitions, causing a faster downward cascade to 
the n=2 state. In other words, once an electron is captured at, say, n=70, it can 
cascade down to the n=2 state faster than in the equilibrium case because few 
photons are around to photoexcite it. In addition, the faster downward cascade 
rate is faster than the photoionization rate from the upper state, and one might 
view this as radiative decay stealing some of the depopulation "flux" from pho- 
toionization. Both the faster downward cascade and the lower photoionization 
rate contribute to the faster net recombination rate. 

4.6. Helium 

With larger ionization potentials, Hell and Hel recombined before H. They can 
be included in the same system of equations in the same framework that has been 
described in this article. These are less important for the calculation of the CMB 
anisotropies, so have generally been paid less scrutiny than H recombination. See 
Seager et al. (1999) and Seager et al. (2000) for a detailed description of He 
recombination. 
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5. Summary 

Time dependence in plasma codes is straightforward to implement if one already 
understands the set of static nonequilibrium rate equations. The same equations 
that are used in static nonequilibrium plasma codes can be used in an integration 
scheme to follow the time evolution of the number densities. The temperature 
equation, with heating and cooling processes should also be evolved with time 
if the relevant heating and cooling processes operate on a timescale longer than 
the physical process that motivates time dependence in the first place. An 
implicit method or other method to treat stiff equations — equations with number 
densities that are changing on very different timescales — is usually necessary. 
Numerical algorithms with adaptive stepsize control and that monitor errors 
are helpful. There can be high gain for this relatively straightforward method. 
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